import sys
!pip install --prefix {sys.prefix} hcp_utils
!pip install --prefix {sys.prefix} nilearn
!pip install --prefix {sys.prefix} nibabel
Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com Requirement already satisfied: hcp_utils in c:\users\adam\anaconda3\lib\site-packages (0.1.0) Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com Collecting nilearn Downloading nilearn-0.9.0-py3-none-any.whl (10.1 MB) Requirement already satisfied: scipy>=1.2 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.6.2) Requirement already satisfied: joblib>=0.12 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.0.1) Requirement already satisfied: scikit-learn>=0.21 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (0.24.1) Requirement already satisfied: pandas>=0.24.0 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.2.4) Requirement already satisfied: requests>=2 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (2.25.1) Requirement already satisfied: numpy>=1.16 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (1.20.1) Requirement already satisfied: nibabel>=2.5 in c:\users\adam\anaconda3\lib\site-packages (from nilearn) (3.2.2) Requirement already satisfied: packaging>=14.3 in c:\users\adam\anaconda3\lib\site-packages (from nibabel>=2.5->nilearn) (20.9) Requirement already satisfied: setuptools in c:\users\adam\anaconda3\lib\site-packages (from nibabel>=2.5->nilearn) (52.0.0.post20210125) Requirement already satisfied: pyparsing>=2.0.2 in c:\users\adam\anaconda3\lib\site-packages (from packaging>=14.3->nibabel>=2.5->nilearn) (2.4.7) Requirement already satisfied: pytz>=2017.3 in c:\users\adam\anaconda3\lib\site-packages (from pandas>=0.24.0->nilearn) (2021.3) Requirement already satisfied: python-dateutil>=2.7.3 in c:\users\adam\anaconda3\lib\site-packages (from pandas>=0.24.0->nilearn) (2.8.1) Requirement already satisfied: six>=1.5 in c:\users\adam\anaconda3\lib\site-packages (from python-dateutil>=2.7.3->pandas>=0.24.0->nilearn) (1.15.0) Requirement already satisfied: certifi>=2017.4.17 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (2021.10.8) Requirement already satisfied: idna<3,>=2.5 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (2.10) Requirement already satisfied: chardet<5,>=3.0.2 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (4.0.0) Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\users\adam\anaconda3\lib\site-packages (from requests>=2->nilearn) (1.26.4) Requirement already satisfied: threadpoolctl>=2.0.0 in c:\users\adam\anaconda3\lib\site-packages (from scikit-learn>=0.21->nilearn) (2.1.0) Installing collected packages: nilearn Successfully installed nilearn-0.9.0 Looking in indexes: https://pypi.org/simple, https://pypi.ngc.nvidia.com Requirement already satisfied: nibabel in c:\users\adam\anaconda3\lib\site-packages (3.2.2) Requirement already satisfied: setuptools in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (52.0.0.post20210125) Requirement already satisfied: numpy>=1.14 in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (1.20.1) Requirement already satisfied: packaging>=14.3 in c:\users\adam\anaconda3\lib\site-packages (from nibabel) (20.9) Requirement already satisfied: pyparsing>=2.0.2 in c:\users\adam\anaconda3\lib\site-packages (from packaging>=14.3->nibabel) (2.4.7)
import hcp_utils as hcp
import numpy as np
import nibabel as nb
import nilearn.plotting as plotting
import matplotlib.pyplot as plt
from scipy import signal
import scipy as sci
import random
# load in subject parcel as template cifti
parcelLoc='C:GroCon.dscalar.nii'
parcel=nb.load(parcelLoc)
parcelCort=parcel.dataobj[:,hcp.struct.cortex]
# load in gradient
gradsf='C:hcp.gradients.dscalar.nii'
grads=nb.load(gradsf)
PG=grads.dataobj[0,hcp.struct.cortex]
vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value C:\Users\adam\anaconda3\lib\site-packages\nibabel\nifti1.py:590: UserWarning: Extension size is not a multiple of 16 bytes; Assuming size is correct and hoping for the best warnings.warn( vox offset (=1.7029e+06) not divisible by 16, not SPM compatible; leaving at current value
# initialize a percentile membership vector
percMemb=np.zeros(len(PG.data))
# get percentiles
for p in range(25):
# every 4 percent bin (25 * 5)
percThresh_L=np.percentile(PG.data,(p*4))
percThresh_U=np.percentile(PG.data,((p+1)*4))
# label cifti verts falling in this perc bin
Indices= np.logical_and(PG.data > percThresh_L, PG.data < percThresh_U)
percMemb[Indices]=p
# plot percentiles/4 to ensure it is working
plotting.view_surf(hcp.mesh.inflated_right, hcp.right_cortex_data(percMemb),
bg_map=hcp.mesh.sulc_right)
# load in template time series to replace data
templateTS=nb.load('C:/Users/adam/template-rest_p2mm_masked.dtseries.nii')
templateData=templateTS.get_fdata()
# extract cortex
templateDataCort=templateData[:,hcp.struct.cortex]
# sampling frequency
TR=.8
fs=(1/TR)
# control for psuedo-high-frequency power from discontinuous (motion-masked) frames
# see https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=2):
# parameters matched to real-data pipeline
nyquist = 0.5 * fs
low = lowcut / nyquist
high = highcut / nyquist
b, a = butter(order, [low, high], btype='band')
return b, a
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return y
# Desired cutoff frequencies (in Hz).
lowcut = 0.008
highcut = 0.09
# create filter
butter_bandpass(lowcut, highcut, fs)
# filter data
templateDataCort = butter_bandpass_filter(templateDataCort, lowcut, highcut, fs)
# extract frequency distribution
# get frequency distribution of all vertices
freqs, density = signal.welch(templateDataCort,fs,axis=0)
# average frequency distribution
avgFreqs=density.mean(axis=1)
plt.plot(freqs,avgFreqs)
plt.xlabel('Hz')
plt.xlim(-0.01,0.25)
plt.show()
# extract cumulative distribution, bin widths (https://stackoverflow.com/questions/13476807/probability-density-function-from-histogram-in-python-to-fit-another-histrogram)
cum_counts = np.cumsum(avgFreqs)
bin_widths = np.repeat(freqs[1] - freqs[0],cum_counts.shape[0])
x = cum_counts*bin_widths
y = freqs[0:]
# interpolate density of frequency distribution
inverse_density_function = sci.interpolate.interp1d(x, y)
# use interpolated density to sample realistic frequencies for simulated waves
randFreqs = np.zeros(100)
for i in range(len( randFreqs )):
u = random.uniform( x[0], x[-1] )
randFreqs[i] = inverse_density_function( u )
# plot randomly sampled values
plt.hist(randFreqs,100)
plt.xlabel('Randomly Sampled Values from Distr: Hz')
plt.xlim(-0.01,0.25)
plt.show()
vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value C:\Users\adam\anaconda3\lib\site-packages\nibabel\nifti1.py:590: UserWarning: Extension size is not a multiple of 16 bytes; Assuming size is correct and hoping for the best warnings.warn( vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value
### replace template data with sine waves
# initialize a new matrix using template data
SimData=np.zeros(templateDataCort.shape)
# initialize tracker of last-filled TR
LF_TR=0
# for 100 iterations, create a sine wave of realistic frequency
# in each percentile bin
for i in range(100):
# get realistic frequency
realisticFreq=randFreqs[i]
# according to this frequency, how many TRs would this wave take to cycle?
numTRs=TR/realisticFreq
# create sine wave of this frequency
realisticPeriod=np.linspace(-np.pi,np.pi,int(np.round(numTRs)))
realisticWave=np.sin(realisticPeriod)
# reverse order of sine wave for EZ interpretation
realisticWave=realisticWave[::-1]
# set starting TR for first bin to be 1 TR after where last wave ended
startTR=LF_TR+1
# Now distribute it over all gradient bins to match a 20-second procession
for p in range(25):
# index vertices in this percentile bin
thisBin=np.where(percMemb==p)
# repeat wave to be consistent across all vertices in this bin
# this is like np.repeat but into a 2d array instead of 1d
binWave=np.tile([realisticWave],(len(thisBin[0]),1))
# get ending TR index
endTR=startTR+len(realisticWave)
# break loop if end TR is beyond length of time series
if endTR > templateDataCort.shape[0]:
break
# extract frames to index into
frames=np.arange(startTR,endTR)
# generate indices that will work to modify SimData
i1,i2 = np.ix_(frames,thisBin[0])
# plop realistic wave into next unused TR
SimData[i1,i2]=np.transpose(binWave)
# delay 1 TR from previous bin
startTR=startTR+1
# record last used TR
LF_TR=startTR+len(realisticWave)
# need dummy subcortex to save as cifti
OutputData=np.zeros(templateTS.shape)
OutputData[:,hcp.struct.cortex]=SimData
# create cifti from dataframe and templates
timeAxis=nb.cifti2.SeriesAxis(start=0,size=SimData.shape[0],step=1)
spaceAxis=templateTS.header.get_axis(1)
new_img=nb.Cifti2Image(OutputData,(timeAxis,spaceAxis))
filename='C:SimulatedPropagations.dtseries.nii'
new_img.to_filename(filename)
# plot PSD of simulated data
freqs, density = signal.welch(SimData,fs,axis=0)
# average frequency distribution
avgFreqs=density.mean(axis=1)
plt.plot(freqs,avgFreqs)
plt.xlabel('Hz')
plt.xlim(-0.01,0.25)
plt.show()